﻿using System;
using System.Collections.Generic;
using System.Text;

namespace SmartMathLibrary.SQLServer.Commands
{
    internal static class SpecialFunctionHelper
    {
        /// <summary>
        /// The incomplete Gammafunction T(a,x).
        /// </summary>
        /// <param name="a">The value for evaluation a.</param>
        /// <param name="x">The value for evaluation x.</param>
        /// <param name="gin">Returns the value for the lower incomplete Gammafunction.</param>
        /// <param name="gim">Returns the value for the incomplete Gammafunction.</param>
        /// <param name="gip">Returns the value for the upper incomplete Gammafunction.</param>
        /// <returns>
        /// The value 0 if the computation was succeed otherwise the calculation failed.
        /// </returns>
        public static int IncompleteGammaFunctionHelper(double a, double x, ref double gin, ref double gim,
                                                        ref double gip)
        {
            int k;
            double ga;

            if ((a < 0.0) || (x < 0))
            {
                return 1;
            }

            double xam = -x + a*Math.Log(x);

            if ((xam > 700) || (a > 170.0))
            {
                return 1;
            }

            if (x == 0.0)
            {
                gin = 0.0;
                gim = GammaFunctionHelper(a);
                gip = 0.0;
                return 0;
            }

            if (x <= 1.0 + a)
            {
                double s = 1.0/a;
                double r = s;

                for (k = 1; k <= 60; k++)
                {
                    r *= x/(a + k);
                    s += r;

                    if (Math.Abs(r/s) < 1e-15)
                    {
                        break;
                    }
                }

                gin = Math.Exp(xam)*s;
                ga = GammaFunctionHelper(a);
                gip = gin/ga;
                gim = ga - gin;
            }
            else
            {
                double t0 = 0.0;

                for (k = 60; k >= 1; k--)
                {
                    t0 = (k - a)/(1.0 + k/(x + t0));
                }

                gim = Math.Exp(xam)/(x + t0);
                ga = GammaFunctionHelper(a);
                gin = ga - gim;
                gip = 1.0 - gim/ga;
            }

            return 0;
        }

        /// <summary>
        /// The Gamma function for real numbers.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of the real gamma function at position x.</returns>
        public static double GammaFunctionHelper(double x)
        {
            //Copyright (c) 2003, C. Bond. All rights reserved.
            int i, k, m;
            double ga = 0;
            double gr = 0;
            double r = 0;
            double z = 0;

            double[] g = {
                             1.0,
                             0.5772156649015329,
                             -0.6558780715202538,
                             -0.420026350340952e-1,
                             0.1665386113822915,
                             -0.421977345555443e-1,
                             -0.9621971527877e-2,
                             0.7218943246663e-2,
                             -0.11651675918591e-2,
                             -0.2152416741149e-3,
                             0.1280502823882e-3,
                             -0.201348547807e-4,
                             -0.12504934821e-5,
                             0.1133027232e-5,
                             -0.2056338417e-6,
                             0.6116095e-8,
                             0.50020075e-8,
                             -0.11812746e-8,
                             0.1043427e-9,
                             0.77823e-11,
                             -0.36968e-11,
                             0.51e-12,
                             -0.206e-13,
                             -0.54e-14,
                             0.14e-14
                         };

            if (x > 171.0)
            {
                return 1e308;
            }

            if (x == (int) x)
            {
                if (x > 0.0)
                {
                    ga = 1.0;
                    for (i = 2; i < x; i++)
                    {
                        ga *= i;
                    }
                }
                else
                {
                    ga = 1e308;
                }
            }
            else
            {
                if (Math.Abs(x) > 1.0)
                {
                    z = Math.Abs(x);
                    m = (int) z;
                    r = 1.0;

                    for (k = 1; k <= m; k++)
                    {
                        r *= (z - k);
                    }

                    z -= m;
                }
                else
                {
                    z = x;
                }

                gr = g[24];

                for (k = 23; k >= 0; k--)
                {
                    gr = gr*z + g[k];
                }

                ga = 1.0/(gr*z);

                if (Math.Abs(x) > 1.0)
                {
                    ga *= r;

                    if (x < 0.0)
                    {
                        ga = -Math.PI/(x*ga*Math.Sin(Math.PI*x));
                    }
                }
            }

            return ga;
        }

        /// <summary>
        /// Computes the Bessel function of the first and second kind.
        /// </summary>
        /// <param name="x">The parameter of the Bessel function.</param>
        /// <param name="j0">Bessel function of 1st kind, 0th order.</param>
        /// <param name="j1">Bessel function of 1st kind, 1st order.</param>
        /// <param name="y0">Bessel function of 2nd kind, 0th order.</param>
        /// <param name="y1">Bessel function of 2nd kind, 1st order.</param>
        /// <param name="j0p">Derivative of Bessel function of 1st kind, 0th order.</param>
        /// <param name="j1p">Derivative of Bessel function of 1st kind, 1st order.</param>
        /// <param name="y0p">Derivative of Bessel function of 2nd kind, 0th order.</param>
        /// <param name="y1p">Derivative of Bessel function of 2nd kind, 1st order.</param>
        /// <returns>Returns 0, if the computation of the Bessel function was correct otherwise, false.</returns>
        public static int BesselFunction(double x, ref double j0, ref double j1, ref double y0, ref double y1, ref double j0p,
                                    ref double j1p, ref double y0p, ref double y1p)
        {
            double M_2_PI = 2 / Math.PI;
            double M_PI_2 = Math.PI / 2;
            double M_PI_4 = Math.PI / 4;
            double eps = 1e-15;
            double el = 0.5772156649015329;
            int k;
            double[] a = {
                             -7.03125e-2,
                             0.112152099609375,
                             -0.5725014209747314,
                             6.074042001273483,
                             -1.100171402692467e2,
                             3.038090510922384e3,
                             -1.188384262567832e5,
                             6.252951493434797e6,
                             -4.259392165047669e8,
                             3.646840080706556e10,
                             -3.833534661393944e12,
                             4.854014686852901e14,
                             -7.286857349377656e16,
                             1.279721941975975e19
                         };

            double[] b = {
                             7.32421875e-2,
                             -0.2271080017089844,
                             1.727727502584457,
                             -2.438052969955606e1,
                             5.513358961220206e2,
                             -1.825775547429318e4,
                             8.328593040162893e5,
                             -5.006958953198893e7,
                             3.836255180230433e9,
                             -3.649010818849833e11,
                             4.218971570284096e13,
                             -5.827244631566907e15,
                             9.476288099260110e17,
                             -1.792162323051699e20
                         };

            double[] a1 = {
                              0.1171875,
                              -0.1441955566406250,
                              0.6765925884246826,
                              -6.883914268109947,
                              1.215978918765359e2,
                              -3.302272294480852e3,
                              1.276412726461746e5,
                              -6.656367718817688e6,
                              4.502786003050393e8,
                              -3.833857520742790e10,
                              4.011838599133198e12,
                              -5.060568503314727e14,
                              7.572616461117958e16,
                              -1.326257285320556e19
                          };

            double[] b1 = {
                              -0.1025390625,
                              0.2775764465332031,
                              -1.993531733751297,
                              2.724882731126854e1,
                              -6.038440767050702e2,
                              1.971837591223663e4,
                              -8.902978767070678e5,
                              5.310411010968522e7,
                              -4.043620325107754e9,
                              3.827011346598605e11,
                              -4.406481417852278e13,
                              6.065091351222699e15,
                              -9.833883876590679e17,
                              1.855045211579828e20
                          };

            if (x < 0.0)
            {
                return 1;
            }

            if (x == 0.0)
            {
                j0 = 1.0;
                j1 = 0.0;
                y0 = -1e308;
                y1 = -1e308;
                j0p = 0.0;
                j1p = 0.5;
                y0p = 1e308;
                y1p = 1e308;

                return 0;
            }

            double x2 = x * x;

            if (x <= 12.0)
            {
                j0 = 1.0;

                double r = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    r *= -0.25 * x2 / (k * k);
                    j0 += r;

                    if (Math.Abs(r) < Math.Abs(j0) * 1e-15)
                    {
                        break;
                    }
                }

                j1 = 1.0;
                r = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    r *= -0.25 * x2 / (k * (k + 1));
                    j1 += r;

                    if (Math.Abs(r) < Math.Abs(j1) * 1e-15)
                    {
                        break;
                    }
                }

                j1 *= 0.5 * x;

                double ec = Math.Log(0.5 * x) + el;
                double cs0 = 0.0;
                double w0 = 0.0;
                double r0 = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    w0 += 1.0 / k;
                    r0 *= -0.25 * x2 / (k * k);
                    r = r0 * w0;
                    cs0 += r;

                    if (Math.Abs(r) < Math.Abs(cs0) * 1e-15)
                    {
                        break;
                    }
                }

                y0 = M_2_PI * (ec * j0 - cs0);

                double cs1 = 1.0;
                double w1 = 0.0;
                double r1 = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    w1 += 1.0 / k;
                    r1 *= -0.25 * x2 / (k * (k + 1));
                    r = r1 * (2.0 * w1 + 1.0 / (k + 1));
                    cs1 += r;

                    if (Math.Abs(r) < Math.Abs(cs1) * 1e-15)
                    {
                        break;
                    }
                }

                y1 = M_2_PI * (ec * j1 - 1.0 / x - 0.25 * x * cs1);
            }
            else
            {
                int kz;

                if (x >= 50.0)
                {
                    kz = 8;
                }
                else if (x >= 35.0)
                {
                    kz = 10;
                }
                else
                {
                    kz = 12;
                }

                double t1 = x - M_PI_4;
                double p0 = 1.0;
                double q0 = -0.125 / x;

                for (k = 0; k < kz; k++)
                {
                    p0 += a[k] * Math.Pow(x, -2 * k - 2);
                    q0 += b[k] * Math.Pow(x, -2 * k - 3);
                }

                double cu = Math.Sqrt(M_2_PI / x);

                j0 = cu * (p0 * Math.Cos(t1) - q0 * Math.Sin(t1));
                y0 = cu * (p0 * Math.Sin(t1) + q0 * Math.Cos(t1));

                double t2 = x - 0.75 * Math.PI;
                double p1 = 1.0;
                double q1 = 0.375 / x;

                for (k = 0; k < kz; k++)
                {
                    p1 += a1[k] * Math.Pow(x, -2 * k - 2);
                    q1 += b1[k] * Math.Pow(x, -2 * k - 3);
                }

                j1 = cu * (p1 * Math.Cos(t2) - q1 * Math.Sin(t2));
                y1 = cu * (p1 * Math.Sin(t2) + q1 * Math.Cos(t2));
            }

            j0p = -j1;
            j1p = j0 - j1 / x;
            y0p = -y1;
            y1p = y0 - y1 / x;

            return 0;
        }

        /// <summary>
        /// Computes the Bessel function of the first and second kind.
        /// </summary>
        /// <param name="n">The order of the Bessel function.</param>
        /// <param name="x">Argument of Bessel function of 1st and 2nd kind.</param>
        /// <param name="jn">Bessel function of 1st kind, orders from 0 to nm.</param>
        /// <param name="yn">Bessel function of 2nd kind, orders from 0 to nm.</param>
        /// <param name="jnp">Derivative of Bessel function of 1st kind, orders from 0 to nm.</param>
        /// <param name="ynp">Derivative of Bessel function of 2nd kind, orders from 0 to nm.</param>
        /// <returns>Returns 0, if the computation of the Bessel function was correct otherwise, false.</returns>
        public static int BesselFunction(int n, double x, double[] jn, double[] yn, double[] jnp, double[] ynp)
        {
            double f = 1;
            double f1;
            int k;
            int nm;
            nm = n;

            if ((x < 0.0) || (n < 0))
            {
                return 1;
            }

            if (x < 1e-15)
            {
                int i;

                for (i = 0; i <= n; i++)
                {
                    jn[i] = 0.0;
                    yn[i] = -1e308;
                    jnp[i] = 0.0;
                    ynp[i] = 1e308;
                }

                jn[0] = 1.0;
                jnp[1] = 0.5;

                return 0;
            }

            BesselFunction(x, ref jn[0], ref jn[1], ref yn[0], ref yn[1], ref jnp[0], ref jnp[1], ref ynp[0], ref ynp[1]);

            if (n < 2)
            {
                return 0;
            }

            double bj0 = jn[0];
            double bj1 = jn[1];

            if (n < (int)0.9 * x)
            {
                for (k = 2; k <= n; k++)
                {
                    jn[k] = 2.0 * (k - 1.0) * bj1 / x - bj0;
                    bj0 = bj1;
                    bj1 = jn[k];
                }
            }
            else
            {
                int m = BesselFunctionHelper1(x, 200);

                if (m < n)
                {
                    nm = m;
                }
                else
                {
                    m = BesselFunctionHelper2(x, n, 15);
                }

                double f2 = 0.0;

                f1 = 1.0e-100;

                for (k = m; k >= 0; k--)
                {
                    f = 2.0 * (k + 1.0) / x * f1 - f2;

                    if (k <= nm)
                    {
                        jn[k] = f;
                    }

                    f2 = f1;
                    f1 = f;
                }

                double cs;

                if (Math.Abs(bj0) > Math.Abs(bj1))
                {
                    cs = bj0 / f;
                }
                else
                {
                    cs = bj1 / f2;
                }
                for (k = 0; k <= nm; k++)
                {
                    jn[k] *= cs;
                }
            }

            for (k = 2; k <= nm; k++)
            {
                jnp[k] = jn[k - 1] - k * jn[k] / x;
            }

            double f0 = yn[0];

            f1 = yn[1];

            for (k = 2; k <= nm; k++)
            {
                f = 2.0 * (k - 1.0) * f1 / x - f0;
                yn[k] = f;
                f0 = f1;
                f1 = f;
            }

            for (k = 2; k <= nm; k++)
            {
                ynp[k] = yn[k - 1] - k * yn[k] / x;
            }

            return 0;
        }

        internal static int BesselFunctionHelper1(double x, int mp)
        {
            int i;
            int nn = 0;
            double a0 = Math.Abs(x);
            int n0 = (int)(1.1 * a0) + 1;
            double f0 = 0.5 * Math.Log10(6.28 * n0) - n0 * Math.Log10(1.36 * a0 / n0) - mp;
            int n1 = n0 + 5;
            double f1 = 0.5 * Math.Log10(6.28 * n1) - n1 * Math.Log10(1.36 * a0 / n1) - mp;

            for (i = 0; i < 20; i++)
            {
                nn = (int)(n1 - (n1 - n0) / (1.0 - f0 / f1));

                double f = 0.5 * Math.Log10(6.28 * nn) - nn * Math.Log10(1.36 * a0 / nn) - mp;

                if (Math.Abs(nn - n1) < 1)
                {
                    break;
                }
                n0 = n1;
                f0 = f1;
                n1 = nn;
                f1 = f;
            }

            return nn;
        }

        internal static int BesselFunctionHelper2(double x, int n, int mp)
        {
            double obj;
            int i, n0, nn = 0;
            double a0 = Math.Abs(x);
            double hmp = 0.5 * mp;
            double ejn = 0.5 * Math.Log10(6.28 * n) - n * Math.Log10(1.36 * a0 / n);

            if (ejn <= hmp)
            {
                obj = mp;
                n0 = (int)(1.1 * a0);

                if (n0 < 1)
                {
                    n0 = 1;
                }
            }
            else
            {
                obj = hmp + ejn;
                n0 = n;
            }

            double f0 = 0.5 * Math.Log10(6.28 * n0) - n0 * Math.Log10(1.36 * a0 / n0) - obj;
            int n1 = n0 + 5;
            double f1 = 0.5 * Math.Log10(6.28 * n1) - n1 * Math.Log10(1.36 * a0 / n1) - obj;

            for (i = 0; i < 20; i++)
            {
                nn = (int)(n1 - (n1 - n0) / (1.0 - f0 / f1));

                double f = 0.5 * Math.Log10(6.28 * nn) - nn * Math.Log10(1.36 * a0 / nn) - obj;

                if (Math.Abs(nn - n1) < 1)
                {
                    break;
                }

                n0 = n1;
                f0 = f1;
                n1 = nn;
                f1 = f;
            }

            return nn + 10;
        }
    }
}